Curva de juros: um conto de três fatores

Author

Victor H. V. Souza

Published

November 20, 2023

1 Preparo e exploração dos dados

rm(list=ls())
library(tidyverse)
library(rb3) # Dados da B3 
library(fixedincome) # Manipulação de dados de renda fixa e juros
library(bizdays) # Calendários
library(plotly)

library(xts)

# options(download.file.method="wininet")
# devtools::install_github('thomasp85/gganimate')
library(gganimate)
library(reactable)

paleta_cores = c('#b8acd1', '#e2d5f1', '#c9e8fd', '#526d85', '#4e4466',
                 '#eccfb0', '#e06666', '#6fa8dc', '#0b5394', '#351c75',
                '#4b2328', '#a86800', '#395a57', '#9ab7c4', '#5e627b')

add_title <- function(x,title_str, subtitle_str = ''){
  x |>
  plotly::layout(title = list(text = paste0('<b>',title_str, '</b>','<br>','<sup>',subtitle_str,'</sup>')))
  }

A ideia do notebook é a partir do dos contratos futuros de DI construir a estrutura a termo de juros da economia brasileira. Por fim, vamos aplicar uma Análise de Componentes Principais e encontrar os fatores de risco latentes na curva de juros. Se você já ouviu falar em Nelson-Siegel, você já deve imaginar quais são.

Na análise, serão usados os ajustes dos contratos negociados na B3 — e tem bastante coisa, de boi-gordo a câmbio, de cupom cambial a café. Os dados podem ser obtidos facilmente por meio do pacote rb3. Vou utilizar também o pacote fixedincome. Ambos foram desenvolvidos pelo Wilson Freitas e facilitam bastante a vida.

# Baixandos os dados da B3 (leva um tempinho, lembre-se de manter o cache para não deixar o TI da bolsa maluco)
# ajustesB3 <- futures_mget(first_date = '2006-01-01', last_date =  Sys.Date(), 
#                   do_cache = T)
# saveRDS(ajustesB3,'data/b3_ajustes.rds')

ajustesB3 = readRDS('data/b3_ajustes.rds')

Obtidos os dados, vamos filtrar o data.frame para os contratos de DI1. Note que o fechamento vem em preço, será necessário transformar em taxa. Deixei a conta explicita.

di1_data = ajustesB3 |>
  filter(commodity == 'DI1') |>
  mutate(date_vencimento = rb3::maturity2date(maturity_code),
         date_vencimento_adj = bizdays::following(date_vencimento, 'Brazil/ANBIMA'), # Ajuste pois o vencimento pode cair em um feriado
         du_ate_vcto = bizdays(refdate, date_vencimento_adj, 'Brazil/ANBIMA'), # Calculando o número de dias úteis até o vencimento do contrato
         tx = ((100000/price)^(1/(du_ate_vcto/252))) -1) |> 
  select(refdate, date_vencimento, date_vencimento_adj, du_ate_vcto, price, tx) |>
   filter(du_ate_vcto > 0)

glimpse(di1_data)
Rows: 150,498
Columns: 6
$ refdate             <date> 2006-01-02, 2006-01-02, 2006-01-02, 2006-01-02, 2…
$ date_vencimento     <date> 2006-02-01, 2006-03-01, 2006-04-01, 2006-05-01, 2…
$ date_vencimento_adj <date> 2006-02-01, 2006-03-01, 2006-04-03, 2006-05-02, 2…
$ du_ate_vcto         <dbl> 22, 40, 63, 81, 124, 188, 249, 311, 373, 437, 499,…
$ price               <dbl> 98587.33, 97464.78, 96079.17, 95023.89, 92627.65, …
$ tx                  <dbl> 0.1769999, 0.1755996, 0.1734998, 0.1720998, 0.1684…

Note que esta é uma estrutura dinâmica. O tempo passa, os contratos vencem. Para continuar o exercício, nós precisamos fixar as maturidades para trabalhar com taxas de 1,2,3 anos.

ultima_refdate_plot = di1_data |>
filter(refdate == max(refdate) | refdate %in% c('2023-09-01', '2023-01-03', '2022-10-04', '2020-04-01')) |>
mutate(refdate = as.factor(refdate)) |> 
ggplot() + 
   geom_point(aes(x = date_vencimento, y = tx, colour = refdate)) +
   geom_line(aes(x = date_vencimento, y = tx, colour = refdate)) +
   theme_minimal() +
   scale_y_continuous(labels = scales::percent) +
   scale_colour_manual('', values = paleta_cores) +
   labs(y = 'Taxa', x = 'Data de Vencimento do Contrato')

ggplotly(ultima_refdate_plot) |>
  add_title(title = 'Curva de Juros', subtitle_str = '% a.a.')

Para isso, é necessário interpolar as curvas — não se assuste, você está apenas ligando os vértices para cada data de referência. Assim, será possível pegar du_ate_vcto fixos de cada refdate. O pacote fixedincome facilita o processo. Vamos criar um objeto SpotRateCurve para cada data de referência, adicionar um método de interpolação e salvar apenas os vértices de interesse.

Quanto aos parâmetros: “discrete” é porque o método de capitalização utilizado é o discreto, “business/252” é porque nossa curva de juros é anualizada em 252 dias úteis e “Brazil/ANBIMA” para utilizar o calendário da ANBIMA para definir os dias úteis.

di1_data_l = di1_data |> 
  arrange(refdate, date_vencimento_adj) |>
  split(di1_data$refdate) |>
  purrr::map(function(x){
    curve = fixedincome::spotratecurve(x$tx, x$du_ate_vcto, refdate = unique(x$refdate),
                               'discrete', 'business/252', 'Brazil/ANBIMA')
    interpolation(curve) <- fixedincome::interp_naturalspline()
    curve
  })

curvaDI1 = tibble::tibble(refdate = names(di1_data_l), curvadi1 = di1_data_l) 

# O legal do pacote é que ele tem umas funções de plot plugadas
# p = curvaDI1 |> filter(refdate == max(refdate)) |>
#     pull(curvadi1) %>%
#     pluck(1) |>
#     fixedincome::ggspotratecurveplot(use_interpolation = TRUE)
# ggplotly(p)

# Agora que temos um estrutura a termo interpolada para cada data de referência, conseguimos por exemplo, pegar a taxa de um ano para cada data de referência
dus = c(
        #21, 42, 63, 126, 
        252, 252*2, 252*3, 252*4, 252*5, 252*6, 252*7, 252*8, 252*9, 252*10
        )

di1_constant_maturity = map_dfr(di1_data_l, 
          .f = function(x){
            tx = x[[dus]] |> as.numeric()
            tibble(refdate = x@refdate,
                   maturity = dus,
                   tx = tx)
          }
)
glimpse(di1_constant_maturity)
Rows: 43,150
Columns: 3
$ refdate  <date> 2006-01-02, 2006-01-02, 2006-01-02, 2006-01-02, 2006-01-02, …
$ maturity <dbl> 252, 504, 756, 1008, 1260, 1512, 1764, 2016, 2268, 2520, 252,…
$ tx       <dbl> 0.1637774, 0.1570996, 0.1539637, 0.1511227, 0.1486126, 0.1461…

Para ilustrar o problema que acabamos de resolver, compare o gráfico abaixo com o primeiro gráfico que nós fixemos. Agora o gráfico não desloca mais, temos as taxas de 252 dias (1 ano), 504 dias (2 anos) e assim sucessivamente para cada ponto no tempo.

E essa foi a evolução da nossa curva de juros de 2020 para cá.

p = di1_constant_maturity |> 
filter(refdate == max(refdate) | refdate %in% c('2023-09-01', '2023-01-03', '2022-01-03', '2021-01-04','2020-04-01', '2019-02-11')) |>
mutate(refdate = as.factor(refdate)) |>
    ggplot() +
    geom_point(aes(x = maturity, y = tx, colour = refdate)) +
   geom_line(aes(x = maturity, y = tx, colour = refdate)) +
   theme_minimal() +
   scale_x_continuous(breaks = dus) +
   scale_y_continuous(labels = scales::percent) +
   scale_colour_manual('', values = paleta_cores) +
   labs(y = 'Taxa', x = 'DU')

ggplotly(p) |>
  add_title(title = 'Curva de Juros', subtitle_str = '% a.a.')

E se a gente animar esse gráfico?

p = di1_constant_maturity |>
    group_by(refdate_month = month(refdate), refdate_year = year(refdate)) |> # Pegando apenas o final de cada mês para não ficar muito pesado
    filter(refdate == max(refdate)) |>
    ungroup() |>
    # mutate(refdate = as.factor(format(refdate, '%b-%y'))) |>
    ggplot() +
      geom_point(aes(x = maturity, y = tx, frame = refdate), size = 2, colour = paleta_cores[1], show.legend = F)+
      geom_line(aes(x = maturity, y = tx, frame = refdate), size = 1, colour = paleta_cores[1], show.legend = F) +
      theme_minimal() +
      scale_x_continuous(breaks = dus) +
      scale_y_continuous(labels = scales::percent, n.breaks = 15) +
      # scale_colour_manual('', values = paleta_cores) +
      labs(y = 'Taxa', x = 'DU')

# Você pode criar um plotly com um seletor de datas
ggplotly(p) |>
animation_opts(1000, easing = 'elastic', redraw = FALSE) |>
animation_button(x = 1, xanchor = 'right', y = 0, yanchor = 'bottom') |>
add_title(title = 'Curva de Juros', subtitle_str = '% a.a.')
# Ou se preferir, você pode criar um GIF
p +
  theme_minimal(base_size = 16) +
  transition_time(refdate) +
  shadow_mark(alpha = 0.1) +
  labs(subtitle = 'Ref. Date: {frame_time}')

E mais gráficos.

# Todos os vértices ao longo do tempo
p = di1_constant_maturity |>
    mutate(maturity = as.factor(maturity)) |>
    ggplot() +
    geom_line(aes(x = refdate, y = tx, colour = maturity)) +
   theme_minimal() +
   scale_y_continuous(labels = scales::percent) +
   scale_colour_manual('', values = paleta_cores) +
   labs(y = 'Taxa', x = 'Data')

ggplotly(p) |>
  add_title(title = 'Curva de Juros', subtitle_str = '% a.a. por vértice')
# Superfície com o Plotly
di1_constant_maturity_xts = di1_constant_maturity |>
  pivot_wider(id_cols = refdate, names_from = maturity, values_from = tx)
di1_constant_maturity_xts = xts(di1_constant_maturity_xts |> select(-refdate), order.by = di1_constant_maturity_xts$refdate)

plot_ly(z = ~di1_constant_maturity_xts, y = index(di1_constant_maturity_xts), x = names(di1_constant_maturity_xts),
         colorbar = list(title = "Taxa")) %>%
  add_surface()  %>%
  layout(scene = list(
           legend = list('Taxa'),
           xaxis = list(title = "DU"),
           yaxis = list(title = "Data"),
           zaxis = list(title = "Taxa")
         )) |>
  add_title(title = 'Curva de Juros', subtitle_str = '% a.a., por maturidade e data de referência')

2 Análise de Componentes Principais (PCA)

Agora que temos tudo pronto, vamos seguir com a Análise de Componentes Principais. A ideia do PCA é que ele é capaz de para uma matriz de n x k, encontrar k-1 vetores independentes que contenham a maior parte da variabilidade encontrada nos vetores originais.

Complicando um pouco, o algoritmo busca pelos autovalores e autovetores da matriz de covariância ou correlação das séries, retornando k-1 vetores linearmente independentes.

# Vamos passar o df para wide e renomear os vértices para ficar mais claro
di1_constant_maturity_wide = di1_constant_maturity |>
  pivot_wider(id_cols = refdate, names_from = maturity, values_from = tx)
colnames(di1_constant_maturity_wide) <-  c('refdate',paste0(as.numeric(
  colnames(di1_constant_maturity_wide[-1])
  )/252, 'Y'))

# Conferindo se tem algum NA nos valores
summary(di1_constant_maturity_wide,digits = 4)
    refdate                 1Y                2Y                3Y         
 Min.   :2006-01-02   Min.   :0.02195   Min.   :0.03230   Min.   :0.04211  
 1st Qu.:2010-05-19   1st Qu.:0.07700   1st Qu.:0.08451   1st Qu.:0.09125  
 Median :2014-09-24   Median :0.10985   Median :0.11365   Median :0.11475  
 Mean   :2014-09-29   Mean   :0.10236   Mean   :0.10631   Mean   :0.10916  
 3rd Qu.:2019-02-07   3rd Qu.:0.12637   3rd Qu.:0.12521   3rd Qu.:0.12497  
 Max.   :2023-11-17   Max.   :0.16400   Max.   :0.17442   Max.   :0.17832  
       4Y                5Y                6Y                7Y         
 Min.   :0.04974   Min.   :0.05417   Min.   :0.05887   Min.   :0.06251  
 1st Qu.:0.09533   1st Qu.:0.09790   1st Qu.:0.09982   1st Qu.:0.10100  
 Median :0.11621   Median :0.11674   Median :0.11732   Median :0.11778  
 Mean   :0.11118   Mean   :0.11238   Mean   :0.11337   Mean   :0.11417  
 3rd Qu.:0.12567   3rd Qu.:0.12581   3rd Qu.:0.12586   3rd Qu.:0.12620  
 Max.   :0.18003   Max.   :0.18099   Max.   :0.18142   Max.   :0.18152  
       8Y                9Y               10Y         
 Min.   :0.06454   Min.   :0.06574   Min.   :0.06658  
 1st Qu.:0.10207   1st Qu.:0.10273   1st Qu.:0.10345  
 Median :0.11809   Median :0.11835   Median :0.11845  
 Mean   :0.11483   Mean   :0.11526   Mean   :0.11562  
 3rd Qu.:0.12660   3rd Qu.:0.12668   3rd Qu.:0.12683  
 Max.   :0.18150   Max.   :0.18150   Max.   :0.18149  
# Qual o desvio padrão de cada vértice da estrutura a termo?
apply(di1_constant_maturity_wide[,-1], 2, sd)
        1Y         2Y         3Y         4Y         5Y         6Y         7Y 
0.03292263 0.02923815 0.02638453 0.02448931 0.02324738 0.02233972 0.02163290 
        8Y         9Y        10Y 
0.02105460 0.02063642 0.02031346 

Assim como foi visto no gráfico, aqui podemos ver que a correlação entre os vértices é bem alta. Logo, é bastante provável que com poucos componentes principais consigamos expressar a maior parte da variabilidade das séries.

Para quem já viu um pouco de finanças isso não é surpresa. Dentre as várias teorias que explicam a estrutura a termo, a Teoria das Expectativas diz que as taxas de juros de X anos nada mais é do que uma combinação da taxa de 1 ano e das taxas de 1 ano esperadas. Isso surge de uma relação de não-arbitragem em que um investidor deve ser indiferente entre investir em um título de 2 anos e investir em um título de 1 ano e depois reinvestir por mais 1 ano.

p = di1_constant_maturity_wide |>
    filter(refdate >= '2015-01-01') |>
    select(-refdate) |>
    corrr::correlate(quiet = T) |>
    corrr::rearrange() |>
    corrr::autoplot()
ggplotly(p) |>
  add_title(title = 'Matriz de Correlação', subtitle_str = 'entre os vértices da curva de juros')

Para performar o PCA, temos a opção de tanto usar a matriz de covariância, como usar a matriz de correlação. Para essa aplicação especifica, vamos utilizar a matriz de correlação. Não queremos que o vértice de maior volatilidade domine sobre os demais. Esse problema fica bem mais evidente em outras aplicações em que as séries são bastante diferentes e com componentes idiossincráticos relevantes (e.g. taxas de câmbio).

# Vamos utilizar o pacote factoextra para performar o PCA
# Coloquei o processo numa função para extrairmos apenas o necessário
get_pca_results <- function(wide_data_frame){
  
  # PCA com normalização
  pca <- wide_data_frame |>
    select(-refdate) |>
    prcomp(scale. = T) # Scale = T utilizamos a matriz de correlação
  
  # Correlações
  pca_corr = factoextra::get_pca(pca)$cor |> as.data.frame.matrix() |>
    select(Dim.1:Dim.10) 
  
  # Contribuição de cada PC 
  pca_cotrib = pca |>
    broom::tidy(matrix = 'd') |>
    filter(PC <= 10) |>
    mutate(PC = as.character(PC)) 

  # Utilizando os loadings ou pesos para recuperar os componentes
  PCs = as.matrix(wide_data_frame |> select(-refdate)) %*% (pca$rotation*-1) |>
    as.data.frame.matrix() |>
    mutate(refdate = wide_data_frame$refdate) |>
    select(refdate, everything())
  
  list_results = list('prcomp.res'  = pca,
       'correlation' = pca_corr,
       'contrib' =  pca_cotrib,
       'PCs' = PCs)
  
  return(list_results)
}
di1_pca = get_pca_results(di1_constant_maturity_wide)

Vamos aos resultados. Primeiramente, sim, apenas 1 componente principal explica 97% da variabilidade em todos os 10 vértices que selecionamos da estrutura a termo. Os 3 primeiros componentes explicam praticamente 100%.

Mas o que tem em cada componente? Pois é, essa é a parte que mais me deixou curioso como um eterno aprendiz e iniciante em 999 coisas, embora para outras pessoas seja óbvio.

Vejamos as correlações dos componentes com as séries originais, bem como os loadings ou auto-vetores de cada componente. O loading nada mais é que a transformação linear que você aplica em cada série para obter o componente principal. Na prática, estamos construindo uma combinação linear de vértices. Cada combinação é linearmente independente da outra.

O primeiro componente (PC1 ou Dim.1) tem alta correlação com todos os vértices e ela é bastante parecida em magnitude. De fato, como podemos ver no gráfico com os loadings, ele é uma combinação linear de pesos iguais de todos os vértices. Uma coisa que vou tentar mostrar posteriormente, é que, na prática, é como se estivessemos fazendo uma média de todos os vértices e extraindo um nível médio da curva naquela data de referência.

O segundo componente (PC2) tem correlação negativa com os vértices curtos e positiva com os vértices mais longos. É como se estivéssemos pegando a parte longa da curva e subtraindo da parte curta. Dito de outra forma, o PC2 é um fator que traz a inclinação da curva de juros.

Por fim, o PC3 ficou um pouco menos óbvio. O plot dos loadings tem quase um formato de U invertido, capturando um terceiro fator que é a curvatura da curva de juros.

De fato, essas ideias estão bem documentadas na literatura de finanças. Com esses 3 fatores, deveríamos ser capazes de reproduzir toda a estrutura a termo da taxa de juros.

# Contribuição, Desvio Padrão por Componente principal
reactable(di1_pca$contrib |> filter(PC %in% as.character(c(1:5))),
          defaultColDef = colDef(format = colFormat(digits = 4)),compact = T, pagination = F, fullWidth = F) |>
  reactablefmtr::add_title('Contribuição e Desvio Padrão por Componente Principal')

Contribuição e Desvio Padrão por Componente Principal

# Correlação entre o Componente Principal e o Vértice da Estrutura a Termo
reactable(di1_pca$correlation[,c(1:5)]*-1,
          defaultColDef = colDef(format = colFormat(digits = 2), style = reactablefmtr::color_scales(di1_pca$correlation*-1)), compact = T, pagination = F, fullWidth = F) |>
  reactablefmtr::add_title('Correlação entre os PCs e os Vértice do DI')

Correlação entre os PCs e os Vértice do DI

# Loadings 
f=paste0(1:10, 'Y')
df_loadings = (di1_pca$prcomp.res$rotation*-1) |> as.data.frame()
df_loadings = df_loadings |> mutate(vertice = rownames(df_loadings)) |>
  select(vertice, everything()) |>
  mutate(vertice = factor(vertice, levels = f))

p=df_loadings |>
  pivot_longer(cols = -1, names_to = 'PC', values_to = 'loading') |>
  filter(PC %in% c('PC1', 'PC2', 'PC3')) |>
  ggplot() +
  geom_bar(stat = 'identity', aes(x = vertice, y = loading), fill = paleta_cores[6]) +
  # geom_hline(aes(yintercept = 0), size = .5, linetype = 'dashed') +
  theme_minimal() +
  facet_wrap(~PC) +
  labs(x = 'Vértice', y = '')
ggplotly(p) |>
  add_title(title = 'Loadings dos componentes principais', subtitle_str = '')

Agora, nós vamos comparar: i) o primeiro componente principal (PC1) com uma métrica de nível da curva de juros; ii) o PC2 com uma medida de inclinação; e por fim, iii) o PC3 com uma métrica de curvatura.

Vamos definir a inclinação como sendo a difernça entre os juros de 10Y e os juros de 1Y. A curvatura será dada por 2*Média(4Y, 5Y, 6Y, 7Y)-(1Y+10Y). Já o nível é simplesmente a média de todos os vértices.

#  Para isso vamos calcular os fatores nível, inclinação e curvatura com os dados originais
fatores_curva = di1_constant_maturity_wide |>
  mutate(nivel = rowMeans(di1_constant_maturity_wide[,-1]),
         inclinacao = `10Y`-`1Y`) |>
  rowwise() |>
  mutate(curvatura = 2*mean(c(`4Y`,`5Y`,`6Y`,`7Y`)) - (`1Y`+`10Y`)) |>
  ungroup() |>
  select(refdate, nivel, inclinacao, curvatura)

# Merge com os componentes principais
fatores_vs_pca = merge(di1_pca$PCs, fatores_curva, by = 'refdate')

fatores_vs_pca_long = fatores_vs_pca |>
  pivot_longer(cols = -1)

ay <- list(
  #tickfont = list(color = "red"),
  overlaying = "y",
  side = "right"
)

# PC1 vs. Nível (Média das Vértices)
plot_ly(data = fatores_vs_pca) %>%
  add_lines(x = ~refdate, y = ~nivel, name = "Nível",type = "scatter", mode = "lines",
            line = list(color = paleta_cores[5])) %>%
  add_lines(x = ~refdate, y = ~PC1, name = "PC1", yaxis = "y2",type = "scatter", mode = "lines", line = list(color = paleta_cores[4])) %>%
  layout(
    yaxis2 = ay,
    xaxis = list(title="Date", ticks=fatores_vs_pca$Date),
    yaxis = list(title = '')
  ) |> add_title(title_str = 'PC1 vs Nível da Curva de Juros')
# PC2 vs. Inclinação (10Y - 1Y)
plot_ly(data = fatores_vs_pca) %>%
  add_lines(x = ~refdate, y = ~inclinacao, name = "Inclinação",type = "scatter", mode = "lines", line = list(color = paleta_cores[5])) %>%
  add_lines(x = ~refdate, y = ~PC2, name = "PC2", yaxis = "y2",type = "scatter", mode = "lines",line = list(color = paleta_cores[4])) %>%
  layout(
    yaxis2 = ay,
    xaxis = list(title="Date", ticks=fatores_vs_pca$Date),
    yaxis = list(title = '')
  ) |> add_title(title_str = 'PC2 vs Inclinação da Curva de Juros',
                   subtitle = 'inclinação (eixo esq.), PC2 (eixo dir.)')
# PC3 vs. Curvatura (Média(2Y,3Y,4Y,5Y,6Y) - (1Y+10Y))
plot_ly(data = fatores_vs_pca) %>%
  add_lines(x = ~refdate, y = ~curvatura, name = "Curvatura",type = "scatter", mode = "lines",line = list(color = paleta_cores[5])) %>%
  add_lines(x = ~refdate, y = ~PC3, name = "PC3", yaxis = "y2",type = "scatter", mode = "lines",line = list(color = paleta_cores[4])) %>%
  layout(
    yaxis2 = ay,
    xaxis = list(title="Date", ticks=fatores_vs_pca$Date),
    yaxis = list(title = '')
  ) |> add_title(title_str = 'PC3 vs Curvatura da Curva de Juros',
                 subtitle = 'curvatura (eixo esq.), PC3 (eixo dir.)')

Como pode ser visto, os resultados são bastante satisfatórios para os dois primeiros componentes. No caso da curvatura, o fit já não é tão bom, sobretudo no passado. É possível melhorar fazendo 2*(3Y) - (1Y+10Y) — replicando um pouco mais fielmente os loadings obtidos para o PC3, entretanto imagino que a ideia de “curvatura” perderia um pouco o sentido. Sendo assim, vamos manter o resultado que encontramos. Um possível motivo para o shape dos loadings do terceiro componente não ter sido o esperado diante da literatura seja as discrepâncias de liquidez entre os vértices.

Enfim, essa é uma das intersecções mais interessantes que tive contato entre técnicas de machine learning e finanças.

Enfim, obrigado se leu até aqui e se você encontrou algum erro ou alguma coisa muito estranha, pode me mandar mensagem ;).